classdef GetXe129FID < Table1.NMRGJobs.NMRGJob
    methods
        function obj = GetXe129FID(nmrg)
            obj@Table1.NMRGJobs.NMRGJob(nmrg, 'GetXe129FID');
            obj.defineParameter();
        end
        
        function obj = defineParameter(obj)
            obj.paramParser=inputParser;
            obj.paramParser.addParameter('duration', 60, @(x) x>0);
            obj.paramParser.addParameter('rate',     2e3, @(x) x>0);
            obj.paramParser.addParameter('targetResolution', 1e-3, @(x) x>0.0)
            obj.paramParser.addParameter('guessFrequency',   253.0, @(x) x>0.0);
            obj.paramParser.addParameter('initAmplitude',    0.050, @(x) x>0.0);
            obj.paramParser.addParameter('oversamfac',       4, @(x) ~rem(x,1));
            obj.paramParser.addParameter('description',      'No description', @ischar);
            obj.paramParser.addParameter('isAOMDark' ,          false, @islogical);
            obj.paramParser.addParameter('isShutterDark' ,      false, @islogical);
        end
        
        %%
        function obj = preprocess(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            
            nmrg.setParametricMagnetometer('frequencyZ', nmrg.SystemInfo.WorkingPoint.mainFrequency,...
                                           'isLockTransverseField', true);
            nmrg.setNMR('frequency129', param.guessFrequency, 'frequency131', param.guessFrequency/3.37340);
            nmrg.by.setDriver( nmrg.ziSlave.LockInChannel(1).setOscillator(1) );
            nmrg.by.setAC('frequency', param.guessFrequency, 'amplitude', param.initAmplitude, 'range', 10.0); pause(2);
            if param.isAOMDark
                nmrg.aom.shutdown();pause(0.5);
            end
            if param.isShutterDark
                beep();pause(1);beep();pause(1);beep();
                input('Please block the laser beam manually.');
            end
            
        end
        
        %%
        function obj = dataAcquisition(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            
            measurement_data.start = datetime('now');
            
            nmrg.ziSlave.disableOutAmplitude();pause(0.1);
            sg = MeasurementComponent.Signal.SignalGroup('Xe129FIDGroup');
            sg.subscribe_signal(nmrg.bTrans, ziDemodTime.X, [2, 1, 1]);
            sg.subscribe_signal(nmrg.Xe129,  ziDemodTime.R, [2, 1, 2]);
            measurement_data.fid = sg.getDemodSignal(param.duration, 'rate', param.rate);
            
            measurement_data.finish = datetime('now');
            obj.dataPackage.measurement_data = measurement_data;
            if param.isAOMDark
                nmrg.aom.setVf(); % use AOM default frequency-control voltage (defined in device .xml file).
            end
            if param.isShutterDark
                beep();pause(1);beep();pause(1);beep();
                input('Please remove the beam block.')
            end
        end
        
        %%
        function obj = analyze(obj)
            param = obj.dataPackage.parameters;
            measurementData = obj.dataPackage.measurement_data;
            
            signal = measurementData.fid;
            fidSignal = signal.getNodeData(1);
            fidEnvelope = signal.getNodeData(2);

            %---------------------------------------------
            idx = 1;
            [xData, yData] = prepareCurveData( fidEnvelope.x(idx:end), fidEnvelope.y(idx:end) );
            ft = fittype( 'exp1' );
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.StartPoint = [max(yData) -mean(xData)];
            
            [fitting.result, fitting.gof] = fit( xData, yData, ft, opts );
            fitting.amplitude = fitting.result.a;
            fitting.tau = abs(1.0/fitting.result.b);
            
            %---------------------------------------------
            [P1,fList] = plomb(fidSignal.y,fidSignal.x,[],param.oversamfac,'psd');
            fitting.fft.value = P1;
            fitting.fft.fList = fList;
            idx = (fList > param.guessFrequency-2.0 & fList < param.guessFrequency+2.0);
            [fpeak, fcenter, fwidth] = findpeaks(fitting.fft.value(idx), fitting.fft.fList(idx), 'NPeaks', 1, 'SortStr', 'descend');

            [xData, yData] = prepareCurveData( fList(and(fList<(fcenter+5*fwidth) , fList>(fcenter-5*fwidth))), ...
                                                sqrt(P1(and(fList<(fcenter+5*fwidth) , fList>(fcenter-5*fwidth)))) );
            
            ft = fittype( 'A/sqrt((x-cf)^2+w^2)+bg', 'independent', 'x', 'dependent', 'y' );
            opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
            opts.Display = 'Off';
            opts.StartPoint = [sqrt(fpeak)*fwidth 0 fcenter fwidth];
            opts.Lower = [0 0 fcenter-1 0];
            opts.Upper = [Inf Inf fcenter+1 0.1];
            opts.TolFun = 1e-16;
            opts.TolX = 1e-16;
            
            [fitting.resultF129, fitting.goff129] = fit( xData, yData, ft, opts );
            fitting.f129 = fitting.resultF129.cf;
            fitting.fft.center = fitting.resultF129.cf;
            fitting.fft.width = fitting.resultF129.w;
            fitting.fft.peak = (fitting.resultF129.A/fitting.resultF129.w+fitting.resultF129.bg);
            
            %-----------------------------------
            obj.dataPackage.analysis_result = fitting;
        end
        
        %%
        function obj = postprocess(obj)
            nmrg = obj.boss;
            param = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis = obj.dataPackage.analysis_result;

            %-----------------------------------
            summary.start = datestr(measuredData.start);
            summary.T2 = analysis.tau;
            summary.f129 = analysis.fft.center;
            summary.finish = datestr(measuredData.finish);
            
            obj.dataPackage.summary = summary;
            nmrg.SystemInfo.Xe129FID = summary;
            
            %-----------------------------------
            getSession;
            sess.addData('GetXe129FID', obj.dataPackage, param.description);
        end
        
        %%
        function obj = presentation(obj, viewer)
            param        = obj.dataPackage.parameters;
            measuredData = obj.dataPackage.measurement_data;
            analysis     = obj.dataPackage.analysis_result;
            summary      = obj.dataPackage.summary;
            
            fidSignal = measuredData.fid.getNodeData(1);
            fidEnvelope = measuredData.fid.getNodeData(2);
            
            if nargin == 1
                fig = figure('Name', sprintf('Xe129FID: [%s]', datestr(measuredData.start)), 'WindowStyle', 'Docked'); 
                ax1=subplot(2, 1, 1, 'Parent', fig);
                ax2=subplot(2, 1, 2, 'Parent', fig);
            else
                Table1.NMRGJobs.NMRGJob.updateViewer(viewer, param, summary);
                ax1=subplot(2, 1, 1,  'Parent', viewer.gui.uiHandle.panelPlot);
                ax2=subplot(2, 1, 2,  'Parent', viewer.gui.uiHandle.panelPlot);
            end
           
            idx = 1;
            plot(ax1, fidSignal.x, fidSignal.y, '-',  fidEnvelope.x,  sqrt(2.0)*fidEnvelope.y/100, 'b.', ...
                fidEnvelope.x(idx:end), sqrt(2.0)*analysis.result(fidEnvelope.x(idx:end))/100, 'r-')
            legend(ax1, 'Envelope', 'Exp. Fitting', 'FID Signal', 'Location', 'best');
            xlim(ax1, minmax(fidSignal.x(:)'));
            xlabel(ax1, 'Time (s)'); ylabel(ax1, 'FID Signal (V)');
            title(ax1, sprintf('\\tau = %4.2f s', analysis.tau), 'Interpreter', 'tex');
            grid(ax1, 'on');
            
            fmin = analysis.fft.center-10.0*analysis.fft.width; fmax = analysis.fft.center+10.0*analysis.fft.width;
            peak_xlist = linspace(fmin, fmax, 500);
            plot(ax2, analysis.fft.fList, analysis.fft.value, 'b.-', peak_xlist, (analysis.resultF129(peak_xlist)).^2, 'r-');
            xlabel(ax2, 'Frequency (Hz)'); ylabel(ax2, 'FFT Abs');
            title(ax2, sprintf('f_0= %6.4f Hz', analysis.fft.center), 'Interpreter', 'tex');
            xlim(ax2, [fmin, fmax]); grid(ax2, 'on');
        end
    end
    %%
    methods(Static)
        function res = compare(objList, varargin)
            p = inputParser;
            p.addParameter('isPlot', true, @islogical);
            p.addParameter('viewer', [], @(x)  isa(x, 'Table1.GUI.DataViewer') || isempty(x));
            p.parse(varargin{:});
            
            res = compare@Table1.NMRGJobs.NMRGJob(objList, 'viewer', p.Results.viewer);

            if p.Results.isPlot
                startTime = arrayfun(@(t) t.start, res.summary, 'UniformOutput', false);
                if isempty(p.Results.viewer)
                    f = figure('Name', sprintf('Xe129FID'));
                    ax1 = subplot(2, 1, 1, 'Parent', f);
                    ax2 = subplot(2, 1, 2, 'Parent', f);
                else
                    ax1 = subplot(2, 1, 1, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                    ax2 = subplot(2, 1, 2, 'Parent', p.Results.viewer.gui.uiHandle.panelPlot);
                end
                
%                 fid = arrayfun(@(d) d.fid.getNodeData(1), res.measuredData, 'UniformOutput', false);
                envelope = arrayfun(@(d) d.fid.getNodeData(2), res.measuredData, 'UniformOutput', false);
                fft_fList = arrayfun(@(d) d.fft.fList, res.analysis, 'UniformOutput', false);
                fft_value = arrayfun(@(d) d.fft.value, res.analysis, 'UniformOutput', false);
                
                cla(ax1); hold(ax1, 'on');
                cellfun(@(enve_k) plot(ax1, enve_k.x, enve_k.y/max(enve_k.y), '.-'), envelope);
                hold(ax1, 'off'); grid(ax1, 'on');
                xlabel(ax1, 'Time (s)'); ylabel(ax1, 'FID Signal (V)');
                legend(ax1, startTime, 'Location', 'bestoutside', 'Interpreter', 'tex');
                
                cla(ax2); hold(ax2, 'on');
                cellfun(@(x, y) plot(ax2, x, y, '.-'), fft_fList, fft_value); 
                hold(ax2, 'off'); grid(ax2, 'on');
                xlabel(ax2, 'Frequency (Hz)'); ylabel(ax2, 'FFT Abs');
            end
        end
    end
end

